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SUMMARY 

We consider a problem which arises in the numerical solution of the compressible two-dimensional or 
axisymmctric boundary-layer equations. Numerical methods for the compressible boundary-layer equa- 
tions are facilitated by transformation from the physical ( x,y ) plane to a computational <£/n) P^ne in 
which the evolution of the flow is "slow" in the time-like E, direction. The commonly used Levy-Lees 
transformation results in a computationally well-behaved problem, but it complicates interpretation of 
the solution in physical space. Specifically, the transformation is inherently nonlinear, and the physical 
wall-normal velocity is transformed out of the problem and is not readily recovered. Conventional 
methods extract the wall-normal velocity in physical space from the continuity equation, using finite- 
difference techniques and interpolation procedures. The present speclrally-accurate method extracts the 
wall-normal velocity directly from the transformation itself, without interpolation, leaving the continuity 
equation free as a check on the quality of the solution. The present method for recovering wall-normal 
velocity, when used in conjunction with a highly-accurate spectral collocation method for solving the 
compressible boundary-layer equations, results in a discrete solution which satisfies the continuity equa- 
tion nearly to machine precision. As demonstration of the utility of the method, the boundary layers of 
three prototypical high-speed flows are investigated and compared: the flat plate; the hollow cylinder; 
and the cone. An important implication for classical linear stability theory is also briefly discussed. 


t This work was concluded partially during the aulhor’s period of tenure as a National Research Council Associate at 
NASA Langley Research Center. Completion of the work was accomplished under NASA Contract NASI -18599, Task 
C-3, by the Theoretical Flow Physics Branch, Fluid Mechanics Division, NASA Langley Research Center, with Analyti- 
cal Services and Materials, Inc. 
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1. INTRODUCTION 

In modem aerodynamics, the boundary-layer approximation is an invaluable tool of widespread 
applicability. Although it is still beyond the capability of existing supercomputers to solve the 
compressible Navier-Stokes equations for complete aerodynamic configurations, it is commonplace for 
engineering purposes to patch inviscid outer solutions to the Euler equations with "inner" solutions to 
the boundary -layer equations to obtain realistic lift and drag estimates. A different application, and the 
motivation for this work, lies in the area of stability and transition, for which solutions to the 
boundary-layer equations provide the mean velocity and temperature distributions necessary for linear 
and nonlinear stability analyses. In this latter context, accuracy is quite important, since, in general, the 
stability of wall-bounded flow is extremely sensitive to variations in the mean. 

The boundary-layer equations define an initial-boundary-value-problcm (IBVP) in which the 
sLreamwise spatial coordinate is lime-like. The solution is obtained by slreamwise marching procedures. 
The equations are extraordinarily "stiff', particularly for high-speed compressible flow. Consequently, 
only implicit marching techniques have met with practical success (for example, see the finite- 
difference method of Harris and Blanchard [1], or the spectral collocation method of Pruett and Street! 
[2]). Depending on the geometry of the flow, the time-like derivative may either enhance or undermine 
the diagonal dominance of the Jacobian used in the iteration procedure. To facilitate the numerical 
solution it is customary to transform the boundary-layer equations from physical (jc,y) space to a com- 
putational (£/n) space in which the time-like derivative has "nice" properties. In the ideal situation, a 
similarity solution exists, and the time-like derivative vanishes identically with the proper similarity 
transformation. Similarity solutions exist, however, only for a limited class of flows (e.g.; flow over a 
flat plate in the absence of a streamwise pressure gradient). For non-similar flows, it is desirable that 
the time-like evolution in the transform plane be "slow", and that the time-like derivative contribute to 
diagonal dominance of the Jacobian. 

One transformation which exhibits these traits for a wide class of boundary-layer flows is that 
known commonly as the Levy-Lees transformation [3]*. Although the Levy-Lees transformation results 
in a computationally well-behaved problem, it complicates interpretation of the results in physical 
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space, relative to the more straightforward transformations used for specialized applications by Duck 
{4], and by Pruett and Streelt [2]. First, it is inherently nonlinear, an additional reason why fully- 
implicit methods are necessary. Second, the physical v -velocity is transformed out of the computa- 
tional problem and is not easily recovered. For some applications, this is not of major concern. For 
example, classical linear stability analyses, which invoke the parallel-flow approximation, disregard the 
wall-normal velocity. It is now recognized, however, that the non-parallel effects on the stability of a 
high-speed flow' can be significant [5], and methods arc being adapted and developed to treat non- 
parallelism. Among these are spatial direct numerical simulation (SDNS) [6], multiple scales analyses 
(MS) [7,8], and a recent scheme based on the parabolized stability equations (PSE) [9,10], Each 
requires accurate determination of wall-normal velocity, and the MS and PSE methods require its gra- 
dient as well. The quality (smoothness) of the solution is of particular importance whenever the appli- 
cation requires differentiated velocities. 

Conventional methods exploit finite-difference techniques and obtain the wall-normal velocity 
from the continuity equation (J. E. Harris, private communication). The method presented here, 
designed as a companion to the spectral collocation method for the compressible boundary-layer equa- 
tions (CBLE) developed by Pruett and Street! [2], enjoys two major advantages over conventional 
approaches. First, the wall-normal velocity is computed to spectral accuracy. Second, the wall-normal 
velocity is extracted directly from the coordinate transformation, leaving the continuity equation avail- 
able as a check on the quality of the solution. Using the method of [2] for the CBLE, and the present 
method to extract wall-normal velocity, we obtain a discrete solution which satisfies the continuity 
equation nearly to machine precision. Moreover, we obtain second derivatives of temperature and velo- 
city distributions which are smooth to at least seven decimal places. 

At the heart of the present method lies the non-trivial evaluation of the quantity T| x . In the next 
section the governing equations and non-dimensionalization are discussed, and the Levy-Lees transfor- 
mation is presented. The third section details the numerical method, focusing on two independent 


♦White [3], however, refers to this as the Ilhngworth-Ixvy-Lees-Dorodnitsyn-Probstein-Elliol transformation, 
mentioning also the contribution of Mangier. 
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derivations for T] it both of which lead to complicated expressions. In the fourth section, in which we 
validate the method, both derivations of r\ x arc shown to give virtually identical numerical results. The 
fifth section provides an application of the method whereby the high-speed boundary layers of a flat 
plate, a hollow cylinder, and a sharp cone are compared, with particular attention to their respective 
wall-normal velocity distributions. An important implication regarding the linear stability of the flow 
along a cone is also discussed. Finally, brief concluding remarks are offered in the last section. 

2. GOVERNING EQUATIONS 

We consider the laminar boundary -layer flow along a two-dimensional or axisymmetric body at 
zero angle of incidence. In regions of the flow in which the boundary-layer approximation is valid, the 
flow is governed by [1] 


T“(^P«) + = 0 

OX dy 


(la) 


du „ du dp 13 

p u— + pv — = + — 

dx dy dx r } dy 


, du 
r'a — - 


(lb) 


3 T _ 37 

pu— + pv— 
ox dy 


dp 

U dx + 


1 d 


+ n 

r 

du 

r*Pr dy 

r J K — 

. d K 



(lc) 



(Id) 


y = \Rey ; v = \Re v 


(lc) 


Dimensionless equations (la,b,c) are derived from the compressible Navier-Stokes equations via the 
boundary-layer approximation [3], and describe, respectively, conservation of mass, stream wise momen- 
tum, and energy. Equation (Id) is the equation of state, and Eq. (le) defines a convenient scaling. 
Here, we assume the fluid to be an ideal gas. For j= 0 and j- 1, Eqs. (1) describe, respectively, two- 
dimensional and axisymmetric boundary-layer flow. 

In Eqs. (1), x is the arc length along the body measured from the stagnation point, y is the wall- 
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normal coordinate, and r = r 0 + y cos(J) is the radial coordinate from the axis of revolution, as shown in 
Fig. 1. In general, r 0 = r 0 (jt) and § = <t>(x). Two axisymmetric bodies of particular interest here are 
the cone and the hollow cylinder, the latter of which is the axisymmetric analog of the flat plate. For 
the cone, <|> is the (fixed) cone half-angle, and r 0 ( x) = *sin({). For the hollow cylinder, r 0 is the (fixed) 
radius, and = 0. 

In Eqs. (1), u and v denote dimensionless velocity components in the x and y dimensions, 
respectively. The remaining dimensionless variables 7\ p, ji, and k are, respectively, the temperature, 
density, viscosity, and thermal conductivity. Lengths are normalized by an arbitrary reference length 
L*, and flow quantities are normalized as follows: 


u 





T = 


r 

t; 





( 2 ) 


Throughout this paper, dimensional quantities are denoted by an asterisk. In Eqs. (2), the reference 
velocity u* and the reference density p* are arbitrary, and the reference temperature, viscosity, and ther- 
mal conductivity are defined as 


*2 



Governing equations (1) are closed by assuming |i(7) and k (T) to vary according to the Sutherland law; 
namely, 


|U= K = 


r^o+c,) 

7+Ci 


C,= 


198.6 °R 


(4) 


The dimensionless parameters which arise as a result of the scalings in Eq. (2) are the Mach 
number A/, Reynolds number Re , the Prandtl number TV, and the ratio of specific heats y. Specifically, 


these are defined as 


- 6 - 



Pr^rL 

Mr‘ 


Pr 



; y 



(5) 


where a* = ^yR g T* is the reference speed of sound, R* g is the ideal gas constant, and C* and C* are, 
respectively, the specific heats at constant pressure and constant volume. For ideal gases, C’ p and C v * 
are constant, and R g = C* - C*. 

The Levy-Lees transformation [3], commonly used to facilitate the numerical solution of Eqs. (1), 
has the following form: 


40 ) = Jp c u e \i e rfrdx (6a) 




'l Re p e u e r{ } 





(6b) 


l = 



1 

1 + 


ycoscj) 

ro 


0 = 0 ) 


0=i) 


(6c) 


In Eqs. (6) and subsequent discussion, subscript "c" denotes edge values, which, in general, depend on 
x . With the further transformations 



(7a) 


V = 


24 

P e u t \i e ri J 
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pVr* 


(7b) 


governing Eqs. (1) take the following form in the transform plane [1]: 
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We note that transformed governing equations (8) make use of 


_G_ = 1 

Pe © 


( 9 ) 


which is obtained from the equation of slate (Id) and the boundary-layer approximation = 0. 

dy 

In the experience of the author, transformation (6) is essential for the accurate numerical solution 
of a non-similar boundary-layer flow in which the streamwise velocity profile becomes "thinner" as one 
proceeds downstream; e.g., boundary-layer flow along a sharp cone. For example, when the spectral 
collocation method of Pruett and Streett [2] using the straightforward transformation of Duck [4] is 
applied for the case of A/ e -6.8 flow over a 7° half-angle cone, the Jacobian of the iteration procedure 
becomes very ill-conditioned. For a highly-resolved grid the condition number can be 0(1O 15 ) or 
larger. In contrast, when transformation (6) is incorporated into the method, the condition number for a 
highly-resolved grid is typically 0(1O 8 ). In the former case, the extreme ill-conditioning prevents full 
convergence of the iteration procedure, resulting in an excessively noisy solution. 


3. NUMERICAL METHODS 


From Eqs. (7b) and (9) one obtains the following expression for the wall-normal velocity v in the 
transform (£,q) plane: 


V&Tl) = 


V24 9 u e li e r$ J V 
Vfo r ’ 24 


jf_ an 

p e dx 


( 10 ) 
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Provided the quantities on the right hand side of Eq. (10) arc known, the determination of v is trivial. 
The difficulty with Eq. (10), however, is that a closed- form expression for r| x is not readily available, 
nor is its form simple. Because of the complicated nature of T[ Xi we present two independent deriva- 
tions which we then use for (numerical) cross-validation. In the first method, r[ x is obtained from the 
Jacobians of transformation (6) and its inverse. In the second method, r[ x is shown to be the soludon 
of a linear Fredholm integral equadon whose coefficients are constructed from information available 
from the solution of Eqs. (8). For both methods, the computadonal effort required to evaluate r\ x (and 
subsequently v) is insignificant relative to that required to solve boundary-layer Eqs. (8). 

Before proceeding, we note that the evaluadon of t in the boundary-layer Eqs. (8), and both 
methods for evaluating rj x , make use of the inverse transformauon corresponding to Eqs. (6). Whereas 
the inverse of (6a) is straightforward and is omitted, the inverse of (6b) requires the solution of a qua- 
dratic equation which results in the following numerical ly-stable form: 


where 


y(£»v) = 1 


— Je^T) = y w 0=0) 

w p e Me i 


r 0 (1 + Vl+C ) 


0 = 1 ) 


(Ha) 


cft.n) - 


2cos<f> y 2 d 
r 0 


(lib) 


Note that the planar 0=0) expression must be evaluated to obtain y w prior to evaluadon of the 
axisymmetric 0=1) expression. 


Method 1 


From Eqs. (6) we obtain 
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X L= dZ L= [dx] 1 

dx dx d 4 


5n 

5? 
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(12b) 


(12c) 

dx dy 54 dx 


To derive Eqs. (12) we have made use of the fact that the product of the Jacobians of transformation 
(6) and its inverse must yield the 2x2 identity matrix. Two of the factors necessary to evaluate (12c) 
are easily obtained from Eqs. (6). They are 


± A 

dx 


P £ M e Pe r 0* 


(13a) 


(i3b) - 

dr\ ^iRe p e u e r* 


The remaining factor is complicated. After straightforward though tedious differentiation of Eq. (11) 
with respect to we obtain 
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( 14 ) 


2D 


r 0 Vl+C 54 
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We note that Leibnitz’ Rule has been exploited for differentiation beneath the integral sign in Eq. (11a). 
As in the case of Eq. (11), one must evaluate the planar 0=0) expression prior to the axisymmetric 
(/=1) expression. The wall-normal velocity v follows immediately from Eqs. (14), (12c), and (10). 

Derivatives of the edge and geometry data t(4)e (p« Me >M* > r o>$) arc evaluated by the chain rule as 
follows: 


dx _ dx_,d 4'_i 
dt, ~ dx^dx 


( 15 ) 
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where is given in Eq. (13a). We defer discussion of the procedure for numerically evaluating the 
integral expressions in Eqs. (11) and (14) until after the the presentation of Method 2. 

Method 2 

First, some mathematical preliminaries are in order. In the steps to follow, we exploit the chain 
rule 


I 


36 = 39 rfg, 36 3p 

dx 3 c, dx 3r| 3 x 


and the integral transformation 


y ti 

[h(y)dy =|A( t 1)-^-</ti 


(16) 


(17) 


where h is an arbitrary function of (£,T|) with t, fixed, and where is given in Eq. (13b). 

dri 

Differentiating transformation (6b) with respect to x , we obtain 


*L (x v) = + 

dx y ,y) p e dx u e dx 


+ j 


p dr 0 
r 0 dx 


J ]_ 'Re PeUed d r < 


2\ dx 


dx 1 0 


Is* 


(18) 


By applying Leibnitz’ Rule to differentiate under the integral sign, and by invoking the product rule, we 
obtain 



(19) 


From definition (6c) and some algebraic manipulation, the following expression for results: 

ax 


3 1 
dx 


= J 


d ( |) 1 dr 0 

dx r 0 dx 


tan^^r 1 - + 


d-0 


( 20 ) 


Expanding the first integral expression on the right hand side of Eq. (19) using Eqs. (17) and (20), we 
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obtain the following expression in terms of (in): 


fl di_ = . 

| 9 dx ^ ^ \ iRe p e u e rl 


, dfy 1 dr o 
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dx r 0 dx 
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( 21 ) 


Similarly, the second integral expression on the right hand side of Eq. (19) is derived to be 


] j _ ae d = _ J !_ 

1 0 2 3x y Re p e u e rh 


d%] 130, } 1 96 3n 
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30 


where chain rule (16) has been invoked for — 


The desired integral form is obtained by combining Eqs. (18-22), whereby 


f(W- 

+ j 


1 dpe 1 dli, 
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p e dx U e dx 
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<U t> 1 dr o 


_ J_A 

2i, dx 

llirfn 


- j tan<t> 
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dt] fl 30 


1 1 


(23) 


We note that certain terms containing — — 7 — have exactly canceled in the derivation of Eq. (23). 

To dx 


Expression (23) is a linear Fredholm integral equation of the form 


T\ 

q (^,ri) = jw (^,n )q (^.'n)^ n + p &*i) 


(24a) 


where 


»«.n) - fj- < M » 

<*> 


p (iti) = / (in + g (\)d (in) - e (in) 


(24d) 
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and where 


1 d P‘ . 1 ^ 1 d\ .. ,d* 

<s£x u e dx 21; dx dx 


,tv . . . rfd> 1 dr 0 

g(t,)=j tan<(>— + — 

dx r 0 dx 


d& Tl) ^ |yr/Tl 


/t , 7 i ae , 

e sT' 71 < 24h > 

We relegate the details of the discretization of Eq. (24a) to the Appendix. We note, however, 
that Eq. (24a) is linear and its matrix analog is N th-order; consequently, the computational effort neces- 
sary to extract v is insignificant relative to that of solving Eqs. (8), which require the iterative solution 
of systems of order 3 M . 

In practice, whether by Method 1 or Method 2, evaluation of tj, occurs along contours of con- 
stant For £ fixed, the integrals which appear in both methods (for example; in Eqs. (13) and (24)) 
assume the form 


i(n) 3 jAODdTl (25) 

o 

where h is again an arbitrary function of T|. It remains to describe their spectrally-accurate numerical 
approximation. For this purpose, let the computational domain 0 < q < q iVMV be partitioned into N 
subintervals such that 0 =q 0 <ni<n 2 < ' • ■ <%=TW,u. At the grid points q„ we have 


JL 0 (n=0) 

*(n«)= 2 A *( T 1*) ; Ai(qJ = (26) 

k=0 < " 

J h(T\)dT\ (« = 1 , 2 ,. ..JV) 

If h n =h(T[ n ), and i„ and Ai n are the discrete approximations of /'(q„) and A/(q„), respectively, then the 
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discrete analog of Eq. (26) is the following: 


», = Z A'* ; 5i = Q/T (27) 

*=0 

In Eq. (27) Q is an jVx(jV+ 1) quadrature operator, and, for convenience, the following vector notation 
is adopted: 


1 

s- 
° . 


i 

■ 

h\ 


^2 


r4 



; M 


i 

r • 


1 


(28) 


Conventions similar to Eqs. (28) hold for all other vector quantities. 

Eq. (27) is general in the sense that Q can represent any quadrature rule. For example, if the tra- 
pezoidal quadrature rule is adopted, then Q is bidiagonal and T approximates i (p) to second-order 
accuracy. For the hydrodynamic stability applications which motivate this work, we are interested in 
attaining highly accurate and smooth solutions. Accordingly, we specialize Q to a collocation method 
based on Chebyshev polynomial approximation, for which Q is a dense matrix and for which spectral 
accuracy is attained. The collocation method applied here to the extraction of the wall-normal velocity 
is described in [2], to which die reader is referred for greater detail. For completeness, we summarize 
briefly below. 

Let continuous non-periodic function h(T\) on [-1,1] be approximated by a finite series expansion 
h N (T\) in an orthogonal basis set of Chebyshev polynomials T n (P), namely 


Mti) = Z^7'„(n) 

n=0 


[■< 


T n W) = cosin cos ’f| 


(29a) 


(29b) 


where coefficients h n are termed the "spectrum" of h N (f\). Using the "natural" Gauss-Lobatto set of 
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collocation points 


an 


n* = cosz n ; z„ = — ; (« = 0,1,2,..JV) 


(30) 


we obtain from Eq. (29) the discrete Fourier cosine transform 


h„ =h N <i n„)= 24C0S-""* 
Jfc=G 


TV 


(31) 


The corresponding discrete inverse transform gives the spectrum h n , namely 


where 


J 2 k= 0 or 
= [1 0<£<W 


It is convenient to express Eq. (32) as a matrix-vector operation 


h = P N Jt 


(32a) 


(32b) 


(33) 


where P N is a dense (TV+l)x(jV+l) matrix whose elements are available by inspection from Eqs. (32). 

Interpreted in the light of transform pair (31) and (32), Eqs. (29) define the spectral interpolation 
polynomial, exact, by definition, at the collocation points. Unlike polynomial interpolating series defined 
on equally spaced intervals, series (29) converges uniformly to h(T\) as N-#>o. Moreover, it can be 
shown that, for continuously differentiable functions h, coefficients h n decay to zero faster than any 
finite power of UN as This is termed "spectral convergence". 

Now, to form quadrature operator Q we consider 


J h(j\)dT\ = E n - / h N <Md1\= J Zh k T k (T\)dT\ = £h k j T k (T\)dT\ 
Vi Vi*=° Vi 


( 34 ) 
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In matrix-vector form, Eq. (34) becomes 


= R N h (35) 

where is an jV-vector and Ryv is a dense Nx(N+l) matrix whose elements are obtained from Eq. 
(34) with the help of Eqs. (29) and (30) as follows: 


R kn=-\ Sim cos (kz )dz ( f = 0 \i£% 

z k~ 1 

From Eqs. (33) and (35) the jVx(JV+l) spectral quadrature operator Q is defined as 


(36) 


Q - iV 


(37) 


We note that a continuous mapping ^^(t)) is necessary to take the natural Chebyshev domain 
— 1<T|<1 onto the computational domain . Such a mapping is also useful to redistribute collo- 


cation points, clustering them in regions of high gradients. 


d 

In practice, the metric — L , computed either 
dr\ 


analytically, or numerically to spectral accuracy, is incorporated directly into quadrature operator Q . 
Finally, we note that Q should be computed initially in double-precision arithmetic to avoid catas- 
trophic loss of significance, The double-precision representation is subsequently truncated to single pre- 
cision for numerical quadrature via Eq. (27). Since Q is computed but once, the additional computa- 
tional effort is minimal. 

We close by summarizing briefly the complete algorithm for extracting wall-normal velocity, as 
integrated into the spectral collocation boundary-layer algorithm of Pruett and Streett [2]. Following 
[2], governing equations (8) are solved for discrete marching steps £ 0 <£i<£2< * ‘ , each corresponding 
to a unique stream wise station. Immediately available from the converged boundary-layer solution at 
each fixed \ are the quantities 0, 0^, 0 n , and t which appear in the integrands of Eqs. (11a), (14), and 
(24). These integrals are evaluated numerically to spectral accuracy following the quadrature procedure 
of Eqs. (27-37). In the discrete approximation, the integral expressions are vector quantities. Also 
needed at each discrete ^ are certain scalar quantities: for example; /(£) and g(^) in Eqs. (24) of 
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Method 2, or the coefficient of y w in E q. (14) of Method 1. These are readily evaluated from the 

discrete edge and geometry data x[%{x)]. Currently, we use cubic splines for smooth interpolation of x 

as well as for computation of derivatives of the form dxidx. Following the evaluation of all necessary 

scalar and vector quantities at fixed die discrete approximation of i\ x is obtained either from Eqs. • 

(12-14) in Method 1, or from integral equation (24) in Method 2, Wall-normal velocity is obtained 

subsequently from Eq. (10) and the boundary-layer solution (V/\0). Despite the awkward nature of the 

expressions which comprise Methods 1 and 2, their evaluations are straightforward, are computationally 

efficient, and involve only the numerical machinery already in place for the solution of boundary -layer 

equations (8). 

4. CODE VALIDATION 

Our purpose here is to offer reasonable validation of the present method. To obtain the results of 
this section, we solve Eqs. (8) by the fully-implicit method of Pruett and Streetl [2], which has been 

I 

modified to incorporate the Levy-Lees transformation (6). We exploit second-order backward 1 

differencing in the time-like dimension, although the method allows up to 5th-order differencing. For 
convenience, the marching scheme uses equally-spaced steps in physical space. A further modification _ 

permits the option of either preconditioned Richardson iteration or Newton iteration within each implicit - 

■ 

marching step. Typically, we iterate using the Richardson scheme until the discrete residual is small ‘ 

enough so that one final (compulationally-intensive) Newton iteration achieves ''full" convergence. Fol- 
lowing convergence of the iteration, the wall-normal velocity is extracted by the spectral collocation 
method presented herein. All computations assume that the wall is adiabatic, although the boundary- 
layer code permits fixed temperature wall conditions as well. The results shown below were obtained 
using 102 (jV= 101) collocation points in the wall-normal T| direction, and constant increments of 
Ax* =0.1 ft. in the marching direction. The following transformation has been used to map from the - 

Chebyshev domain T\ 6 [-1,1] to the computational domain q e[0,Tj Wj4 xl, as discussed in the previous 
section. 
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In Eq. (38), a is a free parameter which controls the strength of stretching. The results shown use 
T] MAX =2l and 0=0.7, chosen after some numerical experimentaUon. For a well-resolved grid, the solu- 
tion is not particularly sensitive to the choice of o. However, it is important that i\ MAX be sufficiently 
large that the far-fteld boundary condition does not "pinch" the boundary-layer flow. Our experience is 
that t\ ma x should be chosen so that y‘OW) is 2 to 3 boundary-layer displacement thicknesses from 
the wall. We note that 102 collocation points are far more than the 40 or so necessary to obtain 3-digit 
"engineering" accuracy. In fact, as will be shown, at this resolution velocity and temperature distribu- 
tions at fixed \ are smooth to 13 digits. 

We take the flow parameters and the geometry of the validation case from the high-speed 
(M„=8.0) wind-tunnel experiment of Stetson et al on a sharp cone at zero angle of incidence [11]. 

These are: 

^ = 7° ; M e — 6.8 ; T* e = 128° R ; = 1.43x10 s ff‘ (39) 

where Re x is the unit Reynolds number based on edge conditions. Except quite near the up, the flow 
on a sharp cone exhibits conditions at the boundary-layer edge which are approximately constant. 
Accordingly, we assume that the edge values remain constant in x , and we set reference values equal to 
their respective edge values; e.g., tf-C We note, however, that both methods for extracting v have 
been validated for the fully general case in which both the geometry and edge data vary with x. 

Figure 2 compares the radial velocity obtained by the present method with the results of a para- 
bolized Navier-Stokes (PNS) calculation using the well-tested code of Korte [12]. The PNS code 
exploits finite-difference techniques and is of second-order accuracy in both the marching (axial) 5 and 
cross-stream (radial) r directions. The velocities («>') computed by the PNS code are in cartesian 
coordinates, and for comparison, the results of the boundary-layer calculation are transformed accord- 
ingly. The PNS scheme is fully explicit and the marching step size is the maximum allowed by stabil- 
ity (CFL) considerations. The PNS computation uses 150 grid points in the radial dimension, approxi- 
mately 30 of which lie within the boundary layer. Al this resolution, the PNS calculation is fairly 
severely under-resolved in the radial direction as indicated by the clearly visible grid-scale oscillation. 
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In fact, as can be inferred from Harris [1], approximately 200 points are necessary to fully resolve the 
boundary layer region to engineering (3-digit) accuracy using 2nd-order finite difference methods. 
Unfortunately, it is impractical for existing computer resources to pack significantly more points within 
the boundary layer since, in the fully explicit PNS code, the marching step size decreases as A r^, 
where A r M}N is the minimum grid interval in the radial dimension. Similar computational barriers were 
encountered when attempting to compare the present method against a compressible Navier-Slokes (NS) 
calculation using the finite-volume code of Jacobs [13], for which it was impractical to resolve the flow 
to the degree desired. Despite these difficulties, the radial velocities obtained from the PNS and 
boundary-layer (BL) computations, which are compared at 5=3.28 ft. (1 meter) in Fig, 2, are in reason- 
able agreement. The maximum velocity, the boundary-layer edge location, and the far-field decay of v' 
all agree well. We note also that the shock obliqueness angle derived by the PNS calculation is about 
1 1 , for which the corresponding post-shock Mach number is 6.8, in concurrence with the experimental 
results of Stetson et al and the edge Mach number used for the BL calculation. We further mention 
that the (under-resolved) NS [13] and PNS [12] calculations are in close agreement. 

The present method has been further validated by Kopriva (in unpublished work), who has 
employed the spectral-collocation BL code to check recent adaptations of his spectral multidomain code 
(inviscid) [14] to incorporate the viscous terms of the CNSE. Agreement between Kopriva’s 
moderately-resolved NS computations and the present method is quite good for test cases which include 
Mach 2 flow over a circular cylinder and Mach 2.2 flow over a flat plate. 

The previous checks against existing codes assure us that there are no grievous errors in the 
current method of extracting wall-normal velocity. We turn now to self-consistency checks. One 
motivation for developing independent methods for evaluating t\ x is to provide a check otherwise una- 
vailable. For the same test case as before, with parameter values given in Eqs. (39), Fig. 3 compares 
T[ x at x =2.0 ft. as computed by Methods 1 and 2. These results appear virtually identical in Fig. 3, 
and, in fact, they agree to at least 1 1 significant digits at every gridpoint. 

As mentioned in the introduction, the current method also leaves the continuity equation available 
as a check on the quality of the discrete solution. For later use, it is preferable to expand continuity 
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Eq. (la), whereby we obtain the following expression, which is valid for flat-plate, hollow cylinder, and 
cone geometry: 


d(P«) , d(P*0 , - P« sin<t> ■ pvcosj), = Q (40) 

dx dy r r 

Figure 4 presents the discrete residual of Eq. (40) at **=3.0 ft. along the cone, computed by summing 
the four terms on the left hand side. Again, the parameter values are those given in Eqs. (39). Deriva- 
tives are evaluated in computational space by the appropriate chain rules, thereby avoiding interpolation 
in physical space. The continuity equation is satisfied to approximately 11 orders of magnitude. That 
there remains considerable structure in the wall-normal distribution of the residuals suggests that, with 
further tuning of o and N , one could likely drive the residual at least another order of magnitude 
toward machine zero (1CT 14 ). 

One final measure of the quality of the spectral numerical method is the decay of the spectra 
(refer to Eq. 32a) as shown at station x* =3.0 ft. along the cone in Fig. 5. The decay of the tempera- 
ture, u -velocity, and v -velocity spectra each by at least 13 orders of magnitude implies that the solution 
is smooth to nearly the full 14-digit precision of the (Cray 2) machine. The linear decay rate on the 
logarithmic scale is indicative of "spectral convergence", by which we mean that truncation error decays 
faster than any finite power of UN as N -**>. 

Our interest in an accurate and smooth solution is not just academic. Analyses of stability based 
on parallel linear stability theory result in eigenvalue problems which require first and second partial 
derivatives with respect to y of both the mean u -velocity and the mean temperature. In addition, ana- 
lyses based on non-parallel theory require the mean v -velocity and its derivatives as well. For AT=101, 
each spectrally-accurate numerical differentiation of the boundary-layer solution results in a loss of 
significance of 2-3 digits due to the unavoidable magnification of round-off errors. Yet because of the 
initially high quality of the solution, we obtain first and second derivatives which are smooth to at least 
10 and 7 digits, respectively. 

Finally, we comment briefly on computational efficiency. For N=101, the spectral collocation 
method of Pruett and Streett [2] modified for the Levy-Lees transformation (6) requires slightly less 
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than one CPU-second per marching step on a Cray 2 supercomputer to solve boundary-layer equations 
(8). Of this second, extraction of the wall-normal velocity by Method 2 consumes a small percentage, 
whereas the computational requirements of Method 1 are virtually negligible. It was originally thought 
that the integral equation method (Method 2) might enjoy some advantages in the distribution of error. 
However, from the data for Fig. 3, it appears that any such advantage is insignificant. 


5. APPLICATION 

We turn now to an application of the method to a problem of both longstanding and immediate 
interest which concerns the stability characteristics of the boundary layers of cylinders and cones rela- 
tive to those of the more commonly studied flat plate. While it is beyond the scope of this paper to 
address stability theory per se t some lighL is shed on stability issues by comparing boundary -layer 
"profiles” for the flat plate, the hollow cylinder, and the cone. 

Here we consider the same flow parameters and cone geometry given in Eqs. (39). Presented for 
comparison, the hollow cylinder can be regarded as the axisymmetric equivalent of the flat plate (that 
is, the boundary layer is assumed to have no thickness at the sharp leading edge). The radius of the 
cylinder is taken to be r 0 =0.3684 ft., equivalent to that of the cone at station jc*= 3.0 ft. For reference, 
Fig. 6 compares the growth of the three boundary layers in terms of displacement thickness 5* as 
defined below [3]: 


5 * 






(41) 


As is well known, the self-similar boundary layer on the flat plate grows proportionately to On 
the other hand, the boundary layer along the cylinder is non-similar because the "curvature" K , 
quantified by the ratio increases in the streamwise direction due to the growth of 5*. Pro- 

vided K<sd, as is the case here, the displacement thickness on the cylinder exhibits growth much like 
that of the flat plate. Also as is well known, the boundary layer on the cone asymptotes to a growth 
rate 1/V3 times that of the flat plate, as predicted by the Mangier transformation [15]. Mangler’s theory 
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ignores transverse curvature and is invalid near the sharp lip of the cone where K is large. As noted by 
Malik and Spall [16], the curvature K along the cone decreases in the streamwise direction, in contrast 
to its increase along the cylinder. As they also note, the effect on stability of the disparity in scales is 
dial the boundary layer on the cone supports higher disturbance frequencies. 

Figures 7a-c compare, respectively, the temperature, streamwise velocity, and wall-normal velo- 
city distributions of the cone, cylinder, and flat-plate boundary layers at x'=3.0 ft. When the wall- 
normal coordinate y* is scaled by the appropriate displacement thickness, the temperature and stream- 
wise velocity distributions of the flat plate, hollow cylinder, and cone "collapse" and are virtually coin- 
cident. There is a slight effect due to transverse curvature which tends to produce a slightly "fuller" 
u -velocity profile and a "thinner" temperature profile with increasing K . For the parameter values of 
the test case, the curvature K values at x * = 3 .0 ft. are 0.0, 0.0744, and 0.0440, for the flat plate, 
cylinder, and cone, respectively. Our principle interest lies in comparison of the wall-normal velocities. 
Whereas the v -velocity in the flat-plate boundary layer is constant at the edge, that of the cylinder 
decays like Mr in the far field. Otherwise, they are qualitatively similar. In contrast to the hollow 
cylinder and flat plate, the cone has v -velocity which changes sign, is negative for large y, and has a 
large, nearly constant gradient in the far field. For comparison Fig. 7c also shows the inviscid solution 
for the wall-normal velocity, obtained from the code of Marconi and co-workers [17]. Note that the 
far-field boundary-layer solution and the inviscid solution have similar trends as anticipated, lending 
additional confidence in the method. 

To highlight the influence of the differences in wall-normal velocity which are due to geometry, it 
is useful to examine the individual contributions of each of the four terms on the left hand side of con- 
tinuity equation (40), as shown in Figs. 8. For the flat plate (Fig. 8a), the third and lourth terms vanish, 
and the first and second terms exactly balance. For the the cylinder (Fig. 8b) there is a small non-zero 
contribution from the fourth term, but the overall picture closely resembles that of the flat plate. In 
contrast, for the cone (Fig. 8c), terms 1 and 2 approximately balance in the near-wall region. However, 
the third term contains a large contribution in the far field which is offset approximately by the far-field 


contribution of the second term. 
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Figure 8c has implications for classical linear stability theory, in which the parallel-flow approxi- 
mation plays an important role. In theory, the parallel-flow approximation arises from the recognition 
that the wall-layer mean flow evolves slowly in x relative to the scale of a typical disturbance 
wavelength, in which case certain terms in the linearized disturbance equations are "small" and can rea- 
sonably be neglected. In practice (i.c., in stability codes), however, the parallel-flow approximation is 
typically implemented by requiring of the mean flow that 

v = 0 and = 0 (42) 

OX 

Condition (42) is self-consistent for the boundary-layer flow along a flat plate or a hollow cylinder in 
the sense that a parallel mean flow can simultaneously satisfy continuity Eq. (40) and approximation 
(42). From Fig. 8c, however, we find that Eqs. (40) and (42) are inconsistent for the cone, since the 
term pusin<{)/r does not vanish under approximation (42). Violation of the continuity equation by the 
assumption of parallelism is particularly noxious since it leads to inconsistency between the conserva- 
tive and nonconservative formulations of the linearized disturbance equations [18]. In fact, Pruett, Ng, 
and Erlebacher [18] have shown rigorously that there is no self-consistent parallel-flow approximation 
for the stability of the flow along a cone. If the stability of a conical flow is properly treated only by 
methods which allow non-parallelism, then the wall-normal velocity assumes importance. Conse- 
quently, the care which has been devoted here to securing an accurate and smooth determination of v is 
not ill spent. 


6, CONCLUSIONS 

The fully implicit, spectral collocation method developed by Pruett and Streett [2] for solution of 
the compressible boundary-layer equations has been extended to incorporate the Levy-Lees [3] transfor- 
mation and spectrally-accurate evaluation of the wall-normal velocity. The Levy-Lees transformation is 
essential to the stability of the numerical method for certain classes of compressible boundary-layer 
flows; however, it renders determination of the wall-normal (v) velocity non-trivial. Two methods of 
determining v are presented and are shown to be numerically equivalent. The generalized algorithm is 
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valid for non-similar, two-dimensional or axisymmetric boundary-layer llow with varying edge condi- 
tions. Computation on a highly resolved mesh results in a discrete solution which satisfies the con- 
tinuity equation nearly to machine precision, while requiring only about one second of CPU time per 
marching step on a Cray 2 supercomputer. Because of its generality, and the accuracy and smoothness 
of the discrete solution, the present method is well suited to providing the mean-flow velocity and tem- 
perature distributions needed for analyses of the stability of compressible boundary-layer flows, whether 
by classical linear stability theory or by recent methods which treat non-parallelism of the mean flow. 
Since mean-flow non -parallel ism can significantly affect the stability of high-speed, wall-bounded flows 
[5], the contribution of wall-normal velocity should not be cavalierly disregarded. 
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APPENDIX 

For completeness, we describe briefly a solution procedure for the discrete analog of integral Eq. 
(24a). For fixed Eq. (24a) has the form 

q(r\) = jw(T])<l(r\)dr\+ p(r\) (Al) 

Note, by inspection of Eqs. (24), 0) = 0. From (Al) we obtain 

<7Cn„) - = J w(f])q(T\)dt\ + (p n - p n -i ) (n = 1,23,. .JV) (A2) 

^-i 

FoUowing the notations regarding Eqs. (26-28), the discrete approximation of Eq. (A2) is 


where W is the NxN diagonal "weighting" matrix 


W = 



w N 


and M is the NxN bidiagonal matrix 


(A3) 


(A4) 


1 

-1 1 
-1 1 


M = 


(A5) 


-1 1 


Q, is the NxN matrix formed by eliminating the first column of the Nx(N+ 1) numerical quadrature 
operator Q (refer to Eq. 37), and is the IV -vector whose components are 

A p n = p n -p n -i ( n = 1,2, 3,.. .AO- Simplifying (A3), we obtain ^ as the solution to the Ath-order 
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linear system 


A = A/?] where A = M - Q i W 


(A6) 
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Radial Velocity 
BL Code vs. PNS Code 
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Fig. 2. Radial velocity at s =3.28 ft. (1 meter). 
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Fig. 4. Residual of the continuity equation. 
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Temperature Distribution 



Fig- 7- Temperature a) ^treamwise velocity b) and wall-normal velo 
city c) distributions at x =3.0 ft. 
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Fig. 7. Temperature a) ^treamwise velocity b) and wall-normal velo- 
city c) distributions at x —3.0 ft. 
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Continuity Equation Terms 

Flat Plate 
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Fig. 8. Terms of the continuity equation compared at station x =3.0: 
a) flat plate; b) cylinder; c) cone. 
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Fig. 8. Terms of the continuity equation compared at station jc * =3.0: 
a) flat plate; b) cylinder; c) cone. 
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Fig. 8. Terms of the continuity equation compared at station * =3.0: 
a) flat plate; b) cylinder; c) cone. 
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